Privratsky-20210215 - Resolution choice
Hélène Fradin Kirshner
02/25/2021
rmarkdown::render(‘./2_Resolution_choice/2_1_Resolution_choice.Rmd’)
Changes in myeloid and kidney cells after CLP - Analysis of 2 x 10X scRNA-seq samples from 2 pools of WT mice (3 Sham + 3 CLP): comparison of gene expression in different cell populations
indir <- "./processedData/1_JP_analyses_results/Rerun_HARDAC_20210216/1_QC_filtering_metrics"
outdir <- "./processedData/2_1_Resolution_choice"
dir.create(outdir, recursive = T)library(Seurat)
filtered <- readRDS(paste0(indir, "/15.filtered.398.rds"))
filtered## An object of class Seurat
## 22399 features across 18055 samples within 1 assay
## Active assay: RNA (22399 features, 0 variable features)Normalize each sample individually and selected 2000 most variable genes between samples
library(cowplot)
list <- SplitObject(filtered, split.by = "sample.id")
list <- lapply(X = list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})
anchors <- FindIntegrationAnchors(object.list = list, dims = 1:20)
integrated <- IntegrateData(anchorset = anchors, dims = 1:20)
DefaultAssay(integrated) <- "integrated"
integrated <- ScaleData(integrated, verbose = T)
integrated## An object of class Seurat
## 24399 features across 18055 samples within 2 assays
## Active assay: integrated (2000 features, 2000 variable features)
## 1 other assay present: RNADefaultAssay(integrated) <- "integrated"
integrated <- RunPCA(integrated, npcs = 50, verbose = FALSE)integrated <- JackStraw(integrated, num.replicate = 100, dims = 50,
verbose = T)
saveRDS(integrated, paste0(outdir, "/1.integrated.rds"))
# integrated <-
# readRDS('./processedData/2_IL1R_KO_vs_ctrl/1.integrated.56.rds')integrated <- ScoreJackStraw(integrated, dims = 1:50)
j <- JackStrawPlot(integrated, dims = 1:50)
j
pdf(paste0(outdir, "/2_JackStrawPlot.pdf"), width = 10, height = 8)
j
dev.off()## png
## 2e <- ElbowPlot(integrated, ndims = 50)
e
pdf(paste0(outdir, "/3_ElbowPlot.pdf"))
e
dev.off()## png
## 2integrated <- RunUMAP(integrated, dims = 1:50)
integrated <- FindNeighbors(integrated, dims = 1:50)
# 0.4-1.2
for (i in seq(0, 3, 0.1)) {
integrated <- FindClusters(integrated, resolution = i, verbose = F)
}
head(integrated[[]])## orig.ident nCount_RNA nFeature_RNA percent.mito sample.id
## AAACCCAAGATGGCGT--C0 C0 572 394 10.839161 C0
## AAACCCAAGCAGTCTT--C0 C0 8811 2172 31.415276 C0
## AAACCCAAGCGAGGAG--C0 C0 12890 3208 21.807603 C0
## AAACCCAAGTAGGGTC--C0 C0 737 405 19.945726 C0
## AAACCCAAGTTTGTCG--C0 C0 634 364 21.451104 C0
## AAACCCACACTAACGT--C0 C0 1932 1151 5.124224 C0
## integrated_snn_res.0 seurat_clusters
## AAACCCAAGATGGCGT--C0 0 16
## AAACCCAAGCAGTCTT--C0 0 3
## AAACCCAAGCGAGGAG--C0 0 2
## AAACCCAAGTAGGGTC--C0 0 14
## AAACCCAAGTTTGTCG--C0 0 40
## AAACCCACACTAACGT--C0 0 0
## integrated_snn_res.0.1 integrated_snn_res.0.2
## AAACCCAAGATGGCGT--C0 2 3
## AAACCCAAGCAGTCTT--C0 0 0
## AAACCCAAGCGAGGAG--C0 0 0
## AAACCCAAGTAGGGTC--C0 1 2
## AAACCCAAGTTTGTCG--C0 3 5
## AAACCCACACTAACGT--C0 3 5
## integrated_snn_res.0.3 integrated_snn_res.0.4
## AAACCCAAGATGGCGT--C0 3 5
## AAACCCAAGCAGTCTT--C0 0 0
## AAACCCAAGCGAGGAG--C0 0 0
## AAACCCAAGTAGGGTC--C0 2 2
## AAACCCAAGTTTGTCG--C0 5 4
## AAACCCACACTAACGT--C0 5 4
## integrated_snn_res.0.5 integrated_snn_res.0.6
## AAACCCAAGATGGCGT--C0 6 6
## AAACCCAAGCAGTCTT--C0 1 2
## AAACCCAAGCGAGGAG--C0 3 3
## AAACCCAAGTAGGGTC--C0 2 1
## AAACCCAAGTTTGTCG--C0 5 4
## AAACCCACACTAACGT--C0 5 4
## integrated_snn_res.0.7 integrated_snn_res.0.8
## AAACCCAAGATGGCGT--C0 5 7
## AAACCCAAGCAGTCTT--C0 1 0
## AAACCCAAGCGAGGAG--C0 2 6
## AAACCCAAGTAGGGTC--C0 3 3
## AAACCCAAGTTTGTCG--C0 4 5
## AAACCCACACTAACGT--C0 4 5
## integrated_snn_res.0.9 integrated_snn_res.1
## AAACCCAAGATGGCGT--C0 8 6
## AAACCCAAGCAGTCTT--C0 2 3
## AAACCCAAGCGAGGAG--C0 4 7
## AAACCCAAGTAGGGTC--C0 5 11
## AAACCCAAGTTTGTCG--C0 3 2
## AAACCCACACTAACGT--C0 3 2
## integrated_snn_res.1.1 integrated_snn_res.1.2
## AAACCCAAGATGGCGT--C0 7 7
## AAACCCAAGCAGTCTT--C0 3 2
## AAACCCAAGCGAGGAG--C0 4 3
## AAACCCAAGTAGGGTC--C0 13 13
## AAACCCAAGTTTGTCG--C0 1 24
## AAACCCACACTAACGT--C0 1 4
## integrated_snn_res.1.3 integrated_snn_res.1.4
## AAACCCAAGATGGCGT--C0 5 8
## AAACCCAAGCAGTCTT--C0 0 0
## AAACCCAAGCGAGGAG--C0 6 9
## AAACCCAAGTAGGGTC--C0 3 4
## AAACCCAAGTTTGTCG--C0 26 26
## AAACCCACACTAACGT--C0 4 3
## integrated_snn_res.1.5 integrated_snn_res.1.6
## AAACCCAAGATGGCGT--C0 7 6
## AAACCCAAGCAGTCTT--C0 0 0
## AAACCCAAGCGAGGAG--C0 2 1
## AAACCCAAGTAGGGTC--C0 13 12
## AAACCCAAGTTTGTCG--C0 26 27
## AAACCCACACTAACGT--C0 3 4
## integrated_snn_res.1.7 integrated_snn_res.1.8
## AAACCCAAGATGGCGT--C0 6 16
## AAACCCAAGCAGTCTT--C0 0 2
## AAACCCAAGCGAGGAG--C0 4 3
## AAACCCAAGTAGGGTC--C0 2 1
## AAACCCAAGTTTGTCG--C0 27 28
## AAACCCACACTAACGT--C0 3 4
## integrated_snn_res.1.9 integrated_snn_res.2
## AAACCCAAGATGGCGT--C0 6 13
## AAACCCAAGCAGTCTT--C0 2 1
## AAACCCAAGCGAGGAG--C0 4 3
## AAACCCAAGTAGGGTC--C0 1 8
## AAACCCAAGTTTGTCG--C0 27 28
## AAACCCACACTAACGT--C0 3 4
## integrated_snn_res.2.1 integrated_snn_res.2.2
## AAACCCAAGATGGCGT--C0 14 9
## AAACCCAAGCAGTCTT--C0 0 4
## AAACCCAAGCGAGGAG--C0 2 1
## AAACCCAAGTAGGGTC--C0 8 12
## AAACCCAAGTTTGTCG--C0 27 30
## AAACCCACACTAACGT--C0 3 0
## integrated_snn_res.2.3 integrated_snn_res.2.4
## AAACCCAAGATGGCGT--C0 12 13
## AAACCCAAGCAGTCTT--C0 7 3
## AAACCCAAGCGAGGAG--C0 1 0
## AAACCCAAGTAGGGTC--C0 16 9
## AAACCCAAGTTTGTCG--C0 32 31
## AAACCCACACTAACGT--C0 0 1
## integrated_snn_res.2.5 integrated_snn_res.2.6
## AAACCCAAGATGGCGT--C0 15 15
## AAACCCAAGCAGTCTT--C0 5 4
## AAACCCAAGCGAGGAG--C0 0 0
## AAACCCAAGTAGGGTC--C0 11 12
## AAACCCAAGTTTGTCG--C0 31 32
## AAACCCACACTAACGT--C0 1 1
## integrated_snn_res.2.7 integrated_snn_res.2.8
## AAACCCAAGATGGCGT--C0 16 16
## AAACCCAAGCAGTCTT--C0 18 3
## AAACCCAAGCGAGGAG--C0 0 0
## AAACCCAAGTAGGGTC--C0 11 14
## AAACCCAAGTTTGTCG--C0 36 34
## AAACCCACACTAACGT--C0 1 1
## integrated_snn_res.2.9 integrated_snn_res.3
## AAACCCAAGATGGCGT--C0 13 16
## AAACCCAAGCAGTCTT--C0 11 3
## AAACCCAAGCGAGGAG--C0 0 2
## AAACCCAAGTAGGGTC--C0 10 14
## AAACCCAAGTTTGTCG--C0 17 40
## AAACCCACACTAACGT--C0 23 0for (i in seq(0.2, 0.3, 0.01)) {
integrated <- FindClusters(integrated, resolution = i, verbose = F)
}# install.packages('clustree')
library(clustree)
c <- clustree(integrated, prefix = "integrated_snn_res.")
c
pdf(paste0(outdir, "/4_clustree.pdf"), width = 17, height = 22)
c
dev.off()## png
## 2# install.packages('clustree')
c <- clustree(integrated, prefix = "integrated_snn_res.", node_colour = "Il6",
node_colour_aggr = "mean")
c
pdf(paste0(outdir, "/5_clustree_Il6.pdf"), width = 17, height = 22)
c
dev.off()## png
## 2# install.packages('clustree')
c <- clustree(integrated, prefix = "integrated_snn_res.", node_colour = "Il10",
node_colour_aggr = "mean")
c
pdf(paste0(outdir, "/6_clustree_Il10.pdf"), width = 17, height = 22)
c
dev.off()## png
## 2Cd11c = Itgax?
# install.packages('clustree')
c <- clustree(integrated, prefix = "integrated_snn_res.", node_colour = "Irf8",
node_colour_aggr = "mean")
c
pdf(paste0(outdir, "/7_clustree_Irf8.pdf"), width = 17, height = 22)
c
dev.off()## png
## 2saveRDS(integrated, paste0(outdir, "/8.integrated.rds"))
# integrated <-
# readRDS('./processedData/2_IL1R_KO_vs_ctrl/1.integrated.56.rds')DefaultAssay(integrated) <- "RNA"
c <- clustree(integrated, prefix = "integrated_snn_res.", node_colour = "Itgax",
node_colour_aggr = "mean")
c
pdf(paste0(outdir, "/9_clustree_Itgax_RNA.pdf"), width = 17,
height = 22)
c
dev.off()## png
## 2MHCII expression to separate from monocytes. Maybe use Ciita ? (Salei paper, p. 16 right before Discussion)
c <- clustree(integrated, prefix = "integrated_snn_res.", node_colour = "Ciita",
node_colour_aggr = "mean")
c
pdf(paste0(outdir, "/10_clustree_Ciita_RNA.pdf"), width = 17,
height = 22)
c
dev.off()## png
## 2Macrophage markers: C1qc, Cd81, Cd74, Apoe (see Zimmerman paper) Dendritic cells: Clec9a, Zbtb46, (Salei, Miller, and Brahler papers)
macrophage.markers <- c("C1qc", "Cd81", "Cd74", "Apoe")
macrophage.markers %in% rownames(integrated)## [1] TRUE TRUE TRUE TRUEmacrophage.markers %in% rownames(integrated@assays$RNA)## [1] TRUE TRUE TRUE TRUEdc.markers <- c("Clec9a", "Zbtb46")
dc.markers %in% rownames(integrated)## [1] TRUE TRUEdc.markers %in% rownames(integrated@assays$RNA)## [1] TRUE TRUEfor (marker in c(macrophage.markers, dc.markers)) {
if (marker %in% macrophage.markers) {
cell.type = "marcophage"
} else {
cell.type = "DC"
}
if (marker %in% rownames(integrated)) {
DefaultAssay(integrated) <- "integrated"
c <- clustree(integrated, prefix = "integrated_snn_res.",
node_colour = marker, node_colour_aggr = "mean")
pdf(paste0(outdir, "/11_clustree_", cell.type, "_", marker,
".pdf"), width = 17, height = 22)
print(c)
dev.off()
} else if (marker %in% rownames(integrated@assays$RNA)) {
DefaultAssay(integrated) <- "RNA"
c <- clustree(integrated, prefix = "integrated_snn_res.",
node_colour = marker, node_colour_aggr = "mean")
pdf(paste0(outdir, "/11_clustree_", cell.type, "_", marker,
"_RNA.pdf"), width = 17, height = 22)
print(c)
dev.off()
}
}Using the stability index to assess clusters The stability index from the SC3 package (Kiselev et al. 2017) measures the stability of clusters across resolutions and is automatically calculated when a clustering tree is built.
c <- clustree(integrated, prefix = "integrated_snn_res.", node_colour = "sc3_stability")
pdf(paste0(outdir, "/12_clustree_integrated_prefix_integrated_snn_res_node_colour_sc3_stability.pdf"),
width = 17, height = 22)
c
dev.off()## png
## 2macrophage.markers <- c("C1qc", "Cd81", "Cd74", "Apoe")
macrophage.markers %in% rownames(integrated)## [1] TRUE TRUE TRUE TRUEmacrophage.markers %in% rownames(integrated@assays$RNA)## [1] TRUE TRUE TRUE TRUEdc.markers <- c("Clec9a", "Zbtb46")
dc.markers %in% rownames(integrated)## [1] TRUE TRUEdc.markers %in% rownames(integrated@assays$RNA)## [1] TRUE TRUEMore markers Podocyte -> Nphs1 / Nphs2 Endo, containing endothelial, vascular, and descending loop of Henle -> Nrp1 Proximal convoluted tubes (PCT) - S1 -> Slc5a2 Proximal tubules - S2 Proximal straight tubules (PST) - S3 -> Atp11a Loop of Henle -> Slc12a1 Distal convoluted tubule (CT) -> Slc12a3 Connecting tubule -> Calb1 CD-PC, collecting duct principal cell -> Aqp2 CD-IC, collecting duct intercalated cell CD-Trans, collecting duct transitional cell Fib, fibroblast -> Plac8 / Pdgfrb Macro, macrophage -> C1qa / Cd68
Neutrophils, eosinophils, and basophils are PMNs Neutro, neutrophil -> S100a8
B lymph, lymphocyte -> Cd79a T lymph, lymphocyte -> Cd4 NK, natural killer cell
Macula densa -> Ptgs2 Granular cell of afferent arteriole -> Ren1 Thin ascending limb -> Clcnka Basophil -> Cd69
markers <- c("S100a8", "Nphs1", "Nphs2", "Nrp1", "Slc5a2", "Atp11a",
"Slc12a1", "Slc12a3", "Calb1", "Aqp2", "Plac8", "Pdgfrb",
"C1qa", "Cd68", "Cd79a", "Cd4", "Ptgs2", "Ren1", "Clcnka",
"Cd69")
for (marker in markers) {
if (marker %in% macrophage.markers) {
cell.type = "marcophage"
} else {
cell.type = "DC"
}
if (marker %in% rownames(integrated)) {
DefaultAssay(integrated) <- "integrated"
c <- clustree(integrated, prefix = "integrated_snn_res.",
node_colour = marker, node_colour_aggr = "mean")
pdf(paste0(outdir, "/13_clustree_", cell.type, "_", marker,
".pdf"), width = 17, height = 22)
print(c)
dev.off()
} else if (marker %in% rownames(integrated@assays$RNA)) {
DefaultAssay(integrated) <- "RNA"
c <- clustree(integrated, prefix = "integrated_snn_res.",
node_colour = marker, node_colour_aggr = "mean")
pdf(paste0(outdir, "/13_clustree_", cell.type, "_", marker,
"_RNA.pdf"), width = 17, height = 22)
print(c)
dev.off()
}
}Session Information
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server release 6.8 (Santiago)
##
## Matrix products: default
## BLAS: /gpfs/fs1/data/omicscore/Privratsky-Privratsky-20210215/scripts/conda/envs/privratsky/lib/libblas.so.3.8.0
## LAPACK: /gpfs/fs1/data/omicscore/Privratsky-Privratsky-20210215/scripts/conda/envs/privratsky/lib/liblapack.so.3.8.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] patchwork_1.1.1 clustree_0.4.3 ggraph_2.0.4 ggplot2_3.3.3
## [5] cowplot_1.1.1 SeuratObject_4.0.0 Seurat_4.0.0
##
## loaded via a namespace (and not attached):
## [1] Rtsne_0.15 colorspace_2.0-0 deldir_0.2-9
## [4] ellipsis_0.3.1 ggridges_0.5.3 spatstat.data_2.0-0
## [7] leiden_0.3.7 listenv_0.8.0 farver_2.0.3
## [10] graphlayouts_0.7.1 ggrepel_0.9.1 RSpectra_0.16-0
## [13] codetools_0.2-18 splines_4.0.3 knitr_1.31
## [16] polyclip_1.10-0 jsonlite_1.7.2 ica_1.0-2
## [19] cluster_2.1.1 png_0.1-7 uwot_0.1.10
## [22] ggforce_0.3.2 shiny_1.6.0 sctransform_0.3.2
## [25] compiler_4.0.3 httr_1.4.2 backports_1.2.1
## [28] Matrix_1.3-2 fastmap_1.1.0 lazyeval_0.2.2
## [31] limma_3.46.0 later_1.1.0.1 tweenr_1.0.1
## [34] formatR_1.7 htmltools_0.5.1.1 tools_4.0.3
## [37] igraph_1.2.6 gtable_0.3.0 glue_1.4.2
## [40] RANN_2.6.1 reshape2_1.4.4 dplyr_1.0.4
## [43] Rcpp_1.0.6 spatstat_1.64-1 scattermore_0.7
## [46] vctrs_0.3.6 nlme_3.1-152 gbRd_0.4-11
## [49] lmtest_0.9-38 xfun_0.20 stringr_1.4.0
## [52] rbibutils_2.0 globals_0.14.0 openxlsx_4.2.3
## [55] mime_0.10 miniUI_0.1.1.1 lifecycle_1.0.0
## [58] irlba_2.3.3 goftest_1.2-2 future_1.21.0
## [61] MASS_7.3-53.1 zoo_1.8-8 scales_1.1.1
## [64] tidygraph_1.2.0 promises_1.2.0.1 spatstat.utils_2.0-0
## [67] parallel_4.0.3 RColorBrewer_1.1-2 yaml_2.2.1
## [70] reticulate_1.18 pbapply_1.4-3 gridExtra_2.3
## [73] rpart_4.1-15 stringi_1.5.3 highr_0.8
## [76] checkmate_2.0.0 zip_2.1.1 Rdpack_2.1
## [79] rlang_0.4.10 pkgconfig_2.0.3 matrixStats_0.58.0
## [82] evaluate_0.14 lattice_0.20-41 ROCR_1.0-11
## [85] purrr_0.3.4 tensor_1.5 htmlwidgets_1.5.3
## [88] labeling_0.4.2 tidyselect_1.1.0 parallelly_1.23.0
## [91] RcppAnnoy_0.0.18 plyr_1.8.6 magrittr_2.0.1
## [94] R6_2.5.0 generics_0.1.0 pillar_1.4.7
## [97] withr_2.4.1 mgcv_1.8-33 fitdistrplus_1.1-3
## [100] survival_3.2-7 abind_1.4-5 tibble_3.0.6
## [103] future.apply_1.7.0 crayon_1.4.1 KernSmooth_2.23-18
## [106] plotly_4.9.3 rmarkdown_2.6 viridis_0.5.1
## [109] grid_4.0.3 data.table_1.13.6 metap_1.1
## [112] digest_0.6.27 xtable_1.8-4 tidyr_1.1.2
## [115] httpuv_1.5.5 munsell_0.5.0 viridisLite_0.3.0writeLines(capture.output(sessionInfo()), "./scripts/2_Resolution_choice/2_1_Resolution_choice.sessionInfo.txt")